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We  present  an  approach  to  simulating  quantum  computation  based  on  a  classical  model  that  directly  imitates 
discrete  quantum  systems.  Qubits  are  represented  as  harmonic  functions  in  a  two-dimensional  vector  space. 
Multiplication  of  qubit  representations  of  different  frequencies  results  in  exponential  growth  of  the  state  space 
similar  to  the  tensor-product  composition  of  qubit  spaces  in  quantum  mechanics.  Individual  qubits  remain 
accessible  in  a  composite  system,  which  is  represented  as  a  complex  function  of  a  single  variable,  though 
entanglement  imposes  a  demand  on  resources  that  scales  exponentially  with  the  number  of  entangled  qubits. 
We  carry  out  a  simulation  of  Shor’s  algorithm  and  discuss  a  simpler  implementation  in  this  classical  model. 
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I.  INTRODUCTION 

Quantum  computation  promises  exponential  speedup  over 
classical  computation  for  certain  problems  such  as  period 
finding  and  quantum  simulation.  Traditional  classical  simu¬ 
lation  of  a  composite  quantum  system  requires  updating  each 
of  the  2n  amplitudes  characterizing  the  state  of  N  qubits, 
according  to  a  Hamiltonian  made  up  of  2NX2N  elements. 
The  exponential  growth  of  the  state  space  with  N  imposes  a 
severe  burden  on  resources  for  this  type  of  simulation.  Find¬ 
ing  classes  of  quantum  computations  that  can  be  simulated 
efficiently  is  an  active  field  of  research  [1]. 

In  quantum  mechanics,  the  qubits  that  comprise  a  com¬ 
posite  system  remain  accessible,  and  it  is  through  interac¬ 
tions  with  individual  and  pairs  of  qubits  that  computation  is 
implemented.  For  example,  a  single-qubit  transformation  af¬ 
fects  all  2n  computational  basis  states  of  an  A'-qubit  system. 
It  therefore  seems  that  one  of  the  features  of  quantum  sys¬ 
tems  that  enables  more  efficient  computing  is  the  ability  to 
harness  the  degrees  of  freedom  of  a  2  VD  vector  space  by 
interacting  with  only  N  qubits. 

Here  we  present  a  classical  model  of  discrete  quantum 
systems  that  is  based  on  representing  individual  qubits  and 
transformations  applied  to  them.  This  enables  us  to  address 
specific  qubits  in  a  composite  system,  the  state  of  which  is 
represented  as  a  complex  function  of  a  single  variable,  and 
thereby  directly  replicate  the  steps  in  an  algorithm  as  they 
would  be  implemented  in  a  quantum  system.  While  the  re¬ 
sources  required  to  carry  out  a  computation  exactly  are  com¬ 
parable  to  other  methods,  this  model  may  be  compatible  with 
approximations  that  would  enable  simulating  more  qubits 
than  currently  is  feasible.  Furthermore,  the  approach  of 
building  a  classical  model  based  on  imitating  quantum  sys¬ 
tems  could  offer  an  opportunity  to  gain  insight  into  the  dif¬ 
ference  in  computational  power  of  classical  and  quantum  ar¬ 
chitectures. 


II.  REVIEW  OF  QUANTUM  COMPUTATION 

A  qubit  can  be  realized  with  any  two-state  quantum  sys¬ 
tem  that  can  be  prepared  in  a  general  superposition  of  basis 
states  of  a  two-dimensional  (2D),  complex  vector  space.  For 
multiple  uncoupled  qubits,  the  state  of  the  composite  system 


is  given  by  the  tensor  product  of  the  individual  states,  and 
the  number  of  computational  basis  states  grows  exponen¬ 
tially  with  the  number  of  qubits.  For  example,  the  state  of  an 
A'-qubit  system,  with  each  qubit  in  an  equal  superposition  of 
computational  basis  states  |0)  and  |l),  is  given  by  the  state 
vector, 

-  ;  (|o>  +  |i))1®(|o>  +  |i»2®(|o) 

+  |l))3®  ...  ®(|0)+|l))jy 
1  \N 

(|00...00>+|00...01> 

+  |00...  10)+  ...  +  111 ...  11)).  (1) 

Application  of  a  unitary  operation  to  a  qubit  that  is  part  of  a 
composite  system,  which  constitutes  a  step  in  an  algorithm, 
affects  all  states  in  the  superposition  simultaneously,  illus¬ 
trating  the  massive  parallelism  inherent  in  quantum  compu¬ 
tation. 

Any  quantum  algorithm  can  be  approximated  arbitrarily 
closely  using  just  single  qubit  operations  and  a  generic  two- 
qubit  interaction,  such  as  the  controlled-NOT  (cnot)  gate  [2]. 
The  effects  of  these  operations  can  be  visualized  as  rotations 
and  inversions  of  the  2JV-dimensional  quantum-computer 
state  vector.  The  state  in  Eq.  (1),  constructed  as  a  product  of 
individual  qubit  states,  is  a  special  case.  Almost  all  of  the 
states  the  system  can  occupy  in  its  vector  space  will  be  non- 
separable,  implying  that  entanglement  is  required  for  general 
computation. 

Qubits  and  operations  on  them  are  subject  to  perturbations 
from  the  environment  and  experimental  imperfections.  In 
general,  it  is  believed  that  errors  due  to  decoherence  grow 
exponentially  with  the  number  of  qubits  in  a  system  [3]. 
Realizable  quantum  computation  relies  on  the  ability  to  di¬ 
minish  the  effects  of  these  errors — quantum  error  correction 
and  fault-tolerant  quantum  computation  exploit  entanglement 
and  the  discrete  nature  of  quantum  systems  to  make  this 
possible. 

This  brief  introduction  to  the  fundamental  elements  of 
quantum  computation  emphasizes  the  role  played  by  the 
mathematical  structure  of  the  single-  and  multiple-qubit  vec- 
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tor  spaces.  Quantum  systems  require  this  mathematical  de¬ 
scription,  and  our  classical  model  is  developed  according  to 
this  description  by  building  into  it  the  same  state-space  struc¬ 
ture.  The  result  is  a  method  of  simulation  and  an  architecture 
that  may  offer  insight  into  the  fundamental  advantages  of 
quantum  systems  for  computing. 

III.  INDIVIDUAL  QUBITS 

For  the  representation  of  a  qubit  we  use  a  harmonic  func¬ 
tion  with  frequency  u>.  Orthogonal  functions  s  i  n  ( tot)  and 
cos(ojf  )  serve  as  convenient  basis  states  and  span  a  2D  vector 
space.  We  assign  these  basis  functions  the  role  of  the  com¬ 
putational  basis  states  of  a  qubit, 

|0)  <=>  sin(cof), 

|l)  <=>  cos(wf).  (2) 

Applying  a  general  unitary  transformation  U  to  a  qubit  ip, 

ip(t)  =  a  sin(wf)  +  (3  cos(wf),  (3) 

requires  isolating  the  coefficients  a  and  /?;  each  coefficient 
can  then  be  multiplied  by  the  corresponding  transformed  ba¬ 
sis  function,  yielding  the  transformed  state: 

U\ip{t)\  =  a£/[sin(<«tf)]  +  /3t/[cos(wf)].  (4) 

Orthogonality  of  the  basis  functions  makes  it  straightforward 
to  isolate  the  coefficients  by  taking  the  inner  product  of  the 
corresponding  function  with  >p{t): 

r'l'nto) 

(O  I 

a=—  sin  (u>t')ip(t')dt' , 

77  J  o 

a  f27r/“ 

f3=—\  cos  (wt')ip(t')dt' .  (5) 

ttJo 

A  general  transformation  is  illustrated  in  Fig.  1(a). 

The  modulus  squared  of  the  coefficient,  or  amplitude,  of  a 
basis  function  gives  the  corresponding  “measurement  prob¬ 
ability.”  The  function  representing  a  qubit  can  be  replaced 
with  one  or  the  other  basis  function  according  to  these  prob¬ 
abilities  in  order  to  represent  the  measurement  process.  The 
ability  to  determine  both  measurement  probabilities  for  a 
given  state  eliminates  the  need  to  introduce  and  carry  along 
normalization  coefficients — the  relative  probabilities  can  be 
determined  at  the  time  of  measurement. 

IV.  COMPOSITE  SYSTEMS 
A.  2  VI)  vector  space 

This  model  can  be  extended  to  composite  systems  by  us¬ 
ing  a  different  frequency  for  the  basis  functions  for  each 
two-state  system,  creating  a  new  2D  vector  space  for  each 
qubit.  The  mapping  of  quantum  states  to  functions  for  com¬ 
posite  systems  becomes 

|0)„<=>sin(co„f), 

|l)„  <=>  cos(«„f),  n  =  0, 1,  (6) 


FIG.  1 .  Schematic  showing  a  linear  transformation  on  single  and 
multiqubit  states,  (a)  For  a  single  qubit,  orthogonality  of  the  basis 
functions  enables  isolation  of  the  coefficients,  which  can  then  be 
multiplied  by  the  transformed  basis  functions,  (b)  For  a  composite 
system,  the  generator  G^"\x,t)  enables  the  functional  form  of 
F'-n\t)  to  be  transferred  to  a  different  variable.  This  procedure  is 
analogous  to  addressing  a  qubit  in  a  quantum  system. 

where  n  refers  to  the  nth  qubit  and  N  is  the  total  number  in 
the  system. 

If  N  single -qubit  functions  in  equal  superpositions  are 
multiplied,  the  result  is  a  linear  combination  of  2N  different 
products, 

''P’a'M  =  [sin(tuif)  +  cos(&)1r)][sin(w2f) 

+  cos(co2f)] . . .  [sin(wA,f)  +  cos(tuA,f)] 

=  [sin(co1f)sin(w2f) . . .  sin( 

+  [sin(oi1f)sin(w2f) . . .  cos(c%f)]  +  . . . 

+  [cos(m1f)cos(cu2f) . ..  cos(c%t)].  (7) 

This  is  analogous  to  the  tensor-product  state  for  a  composite 
quantum  system  in  Eq.  (1),  with  the  mapping 

\bib2b3  ...bN)^>  /!1(w1f)/z2(w2f)/!3(w3f)  . . .  hN(wNt ) 

=  HN,j{t)  ■  (8) 

Here,  bn  is  the  binary  value  representing  the  state  of  the  nth 
qubit  of  a  quantum  system,  hn  is  the  basis  function  (sine  or 
cosine)  representing  the  nth  qubit  in  our  model,  and  HNj(t)  is 
the  yth  of  the  2'v  combinations  of  products  of  h„. 

The  functions  II N(  1)  that  naturally  arise  when  representing 
composite  systems  look  like  A-qubit  computational  basis 
states,  and  we  would  like  to  determine  whether  they  too  span 
a  state  space  that  grows  exponentially  with  qubit  number. 
This  will  be  the  case  if  all  of  the  2N  functions  are  orthogonal. 
It  is  easy  to  see  by  expanding  the  products  of  harmonic  func¬ 
tions  in  terms  of  sum  and  difference  frequencies  that,  for 
N  qubits,  the  HN{t)  are  comprised  of  2N~l  Fourier  fre¬ 
quencies  fI/=S^!=1cr/„w,1,  where  07,,  is  1  or  -1.  The  0/ 
are  all  multiples  of  a  fundamental  frequency  given  by 
the  greatest  common  divisor  of  the  qubit  frequencies. 
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ITiund=gc<-Kw]  ,co2, ,  ojv)  [4].  The  orthogonality  of  the  in¬ 
dividual  Fourier  components  can  be  shown  to  lead  to  or¬ 
thogonality  of  the  HN{t)  for  few  qubits;  the  three -qubit  case 
is  demonstrated  in  the  Appendix. 

1.  Orthogonality  of  the  HN(t ) 

To  demonstrate  orthogonality  in  the  general  case  when  the 
2n~x  Fourier  frequencies  are  unique,  we  consider  the  inner 
product  between  two  /V-qubit  functions  /7v>/(  t) 
=  hj  l(wlt) . . .  h]  N(u>Nt)  and  HNk(t)  =  hk  l(wlt) . . .  hkN{wNt), 

r  ^^fund 

HmJJ)  '  HNk(t)  I  [/^,i(wif)^A-,i(wif)] 

do 

'X[hj2{(02t)hk2{a>2t)\  ■  ■  . 

X[hjN(wNt)hkNi(oNt)]dt,  (9) 

where  factors  for  a  given  qubit  have  been  grouped  together 
[5].  Each  product  hj  n(a)nt)hk  n(a>nt)  can  be  written  as  a  func¬ 
tion  of  frequency  2  a>„,  as  either  ^sin(2a>„f)  or 
^[1  ±cos(2«„f)].  The  integrand  in  Eq.  (9)  can  then  be  seen 
to  consist  of  three  types  of  terms.  First,  there  will  always  be 
a  term  that  is  a  product  of  a  function  for  each  qubit, 
/!1(2w,f)...  hN(2a>Nt).  This  can  be  written  as  a  sum  of  Fourier 
components  with  frequencies  that  are  twice  those  of  the 
HN{t)  and  are,  therefore,  also  multiples  of  (lfund.  Because 
integration  of  a  harmonic  function  over  an  integral  multiple 
of  its  period  yields  zero,  this  term’s  contribution  to  the  inner 
product  Hn j{t)-HNk{t)  vanishes. 

Second,  there  can  be  terms  that  only  include 
factors  for  some  qubits,  such  as 
/i1(2ft)1t).../i„_1(2w„_1t)/i„+i(2w„+|t) . .  .hN(2toNt),  which 
arises  when  hjn=hkn.  The  Fourier  frequencies  for  these 
terms  are  also  multiples  of  nfund — the  fundamental  fre¬ 
quency  is  the  greatest  common  divisor  of  the  set  of  all  qubit 
frequencies,  and  necessarily  divides  any  subset  of  them. 
These  terms  therefore  also  vanish  when  integrated  over  an 
interval  of  27r/Hfund. 

Finally,  the  integrand  in  Eq.  (9)  can  have  a  term  of  unity 
(times  1/2^).  This  only  occurs  when  all  qubit  functions  are 
the  same  for  HNj(t)  and  HNk(t ),  i.e.,  HNj(t)  =  HNk(t). 

2.  Redundant  frequencies 

This  argument  is  only  valid  if  all  of  the  Fourier  frequen¬ 
cies  (1/  are  unique.  If  there  are  at  least  two  combinations  of 
qubit  frequencies  that  give  the  same  Fourier  frequency,  we 
can  write  2^=1cr„ ju>n=^=lan  mton,  for  f Terms  that 
enter  this  equation  with  the  same  sign  for  each  Fourier  fre¬ 
quency  cancel,  and  the  remaining  terms  give  <oa+  atb+ . . . 
=  &>„+wg+...,  where  the  frequencies  have  been  arranged  so 
that  all  signs  are  positive.  When  considering  all  inner  prod¬ 
ucts  HNj{i)-HNk(i),  all  possible  combinations  of  harmonic 
factors  in  the  integrand  in  Eq.  (9)  will  arise;  for  some  inner 
product,  there  will  be  a  term  in  the  integrand  like 
h(2u>at)h(2cobt) . .  .h(2(oj)h(2a>pt). ..,  with  Fourier  frequen¬ 
cies  that  include  fl  =  2(u>a+  u>b+  «g-...)  =  0.  For  the 

case(s)  in  which  this  frequency  is  the  argument  of  a  cosine, 
the  constant  term  results  in  a  nonzero  integral,  and  the  dif¬ 
ferent  HN(t)  in  this  case  are  not  orthogonal. 
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Therefore,  for  sets  of  qubit  frequencies  {«„}  that  result  in 
2N~l  unique  Fourier  frequencies,  the  functions  H N{ t)  that 
naturally  arise  when  representing  composite  systems  are  or¬ 
thogonal  and  span  a  2v-dimensional  space.  Unique  Fourier 
frequencies  can  be  ensured  by  using  a  qubit-frequency  defi¬ 
nition  such  as  (on  =  col  2"_1. 

B.  Linear  operations 

In  quantum  computation,  single-qubit  operations  along 
with  a  generic  two-qubit  interaction,  such  as  a  CNOT  gate,  are 
universal.  The  general  approach  to  implementing  operations 
in  our  model  of  composite  quantum  systems  is  the  same  as 
for  a  single  qubit — we  need  to  isolate  the  factor  multiplying 
each  basis  function  [sin(«„f)  or  cos(w„f)]  for  a  particular 
qubit.  In  this  case,  these  factors  will  be  expressions  involving 
other  qubit  basis  functions.  Once  isolated,  they  can  be  mul¬ 
tiplied  by  the  transformed  basis  functions  and  these  recom¬ 
bined  to  generate  the  transformed  composite  function. 

Consider  a  general  1T,v(r),  similar  to  Eq.  (7)  but  with  ar¬ 
bitrary  coefficients  for  the  H N(t) .  We  can  write  T'  v  as  a  sum 
of  two  parts,  one  with  terms  that  include  cos(w„f)  and  one 
with  terms  that  include  sin(«„f): 

/2am  \  /2n-  i  \ 

'JfvM  =  (  2  jcos(w„f)  +  (  2  bkH^N(t)Jsm(a)nt) 

=  F{"\t)cos(a)nt)  +  F^)(f)sin(w„f).  (10) 

The  Hk%(t)  are  products  of  harmonic  functions, 

hi{Mit)h2(to2t)  ■  ■  ■  /z„_1(«„_if)/z)!+1(«„+1f) . . .  hN{wNt), 

where  h  is  cosine  or  sine,  and  ak  and  bk  are  coefficients  for 
the  cos(w;,f)  and  sin(w„f)  terms.  F^\t)  and  F^\t)  are  the 
functions  that  we  need  to  be  able  to  isolate  to  apply  a  linear 
transformation  to  qubit  n;  determining  these  functions  can  be 
considered  to  be  “addressing  qubit  n.” 

For  small  N,  a  procedure  similar  to  the  one  for  a  single 
qubit  can  be  adapted.  Multiplication  of  by  the  relevant 
basis  function  for  the  qubit  leads  to  a  different  spectrum  for 
terms  containing  that  basis  function.  These  different  frequen¬ 
cies  could  be  selected  from  the  terms  containing  the  orthogo¬ 
nal  basis  function,  enabling  the  qubit  to  be  addressed.  As  the 
qubit  number  grows,  however,  the  number  and  density  of 
frequencies  grow  dramatically,  making  this  process  unfea¬ 
sible. 

A  more  general  procedure  can  be  used  which  determines 
exactly  using  the  orthogonality  of  the  FlN(t).  An 
inner  product  can  be  imposed  between  'I'  v  and  a  projector 
that  forces  all  of  the  terms  with  one  basis  function  for  qubit 
n  to  vanish  while  preserving  the  others.  The  construction  of 
the  projector  for  a  given  system  is  straightforward. 

A  linear  combination  of  all  of  the  HN(t)  for  a  system  of  N 
qubits  can  be  generated  by  putting  each  qubit  into  an  equal 
superposition  of  basis  functions  as  in  Eq.  (7).  We  define  a 
similar  function,  the  generator  for  qubit  n,  as  the  product  of 
equal  superpositions  of  computational  basis  states  for  all  qu¬ 
bits  in  the  system  except  qubit  n: 
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G^\t)  =  [sin(o)jf)  +  cos(a)1r)][sin(w2f) 

+  cos(ci>2f)]  ■  ■  ■  [sin(w«-il)  +  COs(<W„_[f)][sin(«„+1f) 

+  cos(«„+1f)] . . .  [sin(wA,f)  +  costw^f)].  (11) 

If  this  is  multiplied  by  cos(oj„f)[sin(oj„f  )].  the  resulting  func¬ 
tion’s  inner  product  with  d'v(f)  gives  the  sum  of  the  ampli¬ 
tudes  of  the  terms  in  but  the  functional  form 

is  lost.  We  can  salvage  the  functional  dependence  by  trans¬ 
ferring  it  to  a  second  variable  introduced  to  exactly  replicate 
the  dependence  on  t : 

=  [sin(w1f)sin(w1x)  +  cos(w1f)cos(<u1x)] 

X[sin(<w2f)sin(<w2x) 

+  cos(co2f)cos(w2x)] . . .  [sin(a>„_1/)sin(<un_1x) 

+  cos(a)„_1f)cos(w,!_1x)][sin(w„+1f)sin(w,i+1x) 

+  cos(co„+1f)cos(<w„+1x)] . . .  [sin(wwf)sin(w;vx) 

+  cos(wwf)cos(o);vx)] 

=  cos(o>i?  -  <u1x)cos(w2f  -  <w2x) . . .  cos(w„_1f 
-  <w„_1x)cos(<w„+1f  -  u>n+\x) . . .  cos (a>Nt  -  u>Nx) . 

(12) 

When  this  generator  is  multiplied  by  cos(<u„r)[sin(w„f)],  we 
get  the  projector  for  F^(t)[F^(t)]: 

P(”\x,t)  =  cos(o )nt)G^\x,t) 


[P("\x,t)  =  sin(w),f)G^,)(x,f)].  (13) 

Taking  the  inner  product  of  and  Pj,")(x,f)[/3^(x,f)], 

integrated  over  t,  gives  us  F<p‘\x)[F<'G  (x)]\ 


ff(x)  =  (2*0^/ 7r) 


Ff\x)  =  (2/vnfund/  7r) 


>1 

>1 


0 

TT/ftf, 


u 

P("\x,t)ii/  N(t)dt 
P^Hxj)^  N(t)dt 


■  (14) 


This  is  an  exact  technique  for  addressing  a  qubit  that  is  part 
of  a  composite  system.  The  rest  of  the  procedure  for  apply¬ 
ing  a  transformation  follows  as  in  the  single-qubit  case  and  is 
illustrated  in  Fig.  1(b). 

Iteration  of  this  technique  of  addressing  qubits  allows  for 
multiple-qubit  gates.  For  example,  for  a  controlled-NOT  gate 
with  qubit  nx  as  the  control  and  n2  as  the  target,  the  gate 
would  begin  with  determination  of  Fy\  which  would  then 
take  on  the  role  of  the  function  'T  for  determining  and 
Fy\  Fy1^  would  be  reconstructed  after  inverting  the  basis 
functions  for  n2,  giving  F^(t)sin(£o2?)+F^(f)cos(w2f).  Fi¬ 
nally,  the  transformed  state  would  be  generated  as 
[^"2)(f)sin((w2f)  +  F'^2)(f)cos(w„f)]cos(wIf)  +  F*"l)(f)sin(<w1f). 
Gates  involving  more  than  two  qubits  can  be  implemented 
by  further  iteration. 

The  measurement  probability  for  a  basis  function  is  deter¬ 
mined  by  fg,slf mdF(n>*(t)-F^(t)dt  [6].  Due  to  the  orthogonal¬ 
ity  of  the  H^t),  all  cross  terms  in  the  inner  product  vanish 
and  the  result  is  the  sum  of  the  moduli  squared  of  the  ampli¬ 


tudes  of  all  of  the  terms  containing  the  corresponding  basis 
function. 


C.  Scaling  of  required  resources 

The  state  of  a  general  composite  system  can  be  repre¬ 
sented  as 


^(wbw2,  ...  ,uN)ip{wNe+\)ip{<0Ne+2)  ■■■  ,  (15) 


where  ip  represents  an  individual,  unentangled  qubit,  and  d' 
characterizes  Ne  entangled  qubits.  While  unentangled  qubits 
can  be  stored  and  processed  individually  and  with  little  over¬ 
head,  the  resources  required  to  exactly  represent  entangled 
qubits  scale  exponentially  with  Ne  [7].  The  qubit  frequencies 
and  the  maximum  Fourier  frequency  can  be  kept  finite,  but 
the  fundamental  frequency  decreases  at  least  exponentially 
with  number.  The  interval  over  which  d'  needs  to  be  defined 
is  given  by  the  integration  interval  required  for  addressing  a 
qubit.  In  Eq.  (14),  the  integration  limit  of  rr/flfund  yields 
exact  values  for  the  F*"-1;  coupled  with  the  necessary  reso¬ 
lution  imposed  by  the  highest  Fourier  frequency,  on  order  of 
flmax/  Qjund  points  are  required  to  define  d'. 

We  can  consider  the  impact  on  addressing  a  qubit — both 
for  a  unitary  transformation  and  for  measurement — of  simply 
truncating  all  of  the  functions  that  arise  in  a  calculation.  For 
unitary  transformations,  we  determine  the  accuracy  of 
F^\r,t)  and  Fy(r,f),  the  functions  determined  by  integrat¬ 
ing  Eq.  (14)  to  t  rather  than  to  7r/(lfund,  by  comparing 
F®(T,f)cos(w„f)+F^(T,f)sin(&i)1f)  to  dr(f).  We  define  a  pa¬ 
rameter  b1",( r)  to  represent  the  error  in 


f 

J  0 


P[n>(x,t)'V(t)dt 


cos(w„x) 


+ 


P{;\x,t)V(t)dt 


sin(w„x) 


2 

dx, 


(16) 


where  d'  is  d'  normalized  over  the  interval  0  to  t,  and  /V(  r) 
is  the  normalization  constant  for  F^(r,  f)cos(«„f) 
+f["\t,  f)sin(<w„f)  over  the  same  interval. 

In  Fig.  2(a)  we  plot  S(  t)  for  the  state  d r(t) 
=  n^1cos(w„r)  +  n^Isin(w„f),  with  a>n=a>l2"~l  and  Ne=5 
through  9.  We  evaluate  S  for  the  case  of  addressing  the  first 
qubit.  The  curves  show  that  determination  of  F^  and  F("\ 
which  is  exact  for  an  integration  limit  of  7r/fllurld,  abruptly 
becomes  less  precise  as  the  integration  interval  is  reduced. 
Any  calculation  involving  many  gates  will  likely  require  a 
value  for  S  on  the  steep  part  of  the  curve,  imposing  an  inte¬ 
gration  interval  that  scales  exponentially  with  Ne. 

To  assess  the  effect  of  truncation  on  measurement  prob¬ 
abilities,  we  truncate  the  integrals  used  to  determine  Fy  and 
Fy  as  above,  and  then  integrate  the  square  of  each.  We  look 
at  the  ratio  of  the  truncated  probability  to  measure  sine  ver¬ 
sus  cosine  as  a  function  of  integration  time: 
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t  (units  of  2n/co) 


2  468  2  468 

10  100 

t  (units  of  27r/co) 


FIG.  2.  (Color  online)  (a)  Semi-logarithmic  plot  of  error  S  versus  integration  time  r  when  addressing  qubit  n=  1  for  the  state  discussed 
in  the  text.  The  equal  spacing  between  the  steep  parts  of  the  curves  indicates  an  exponential  scaling  of  integration  interval  in  order  to  avoid 
significant  errors,  (b)  Plot  of  r  versus  integration  time  t,  for  the  same  state  and  addressing  the  same  qubit.  For  both  plots,  curves  from  left 
to  right  correspond  to  Ne= 5,  6,  7,  8,  and  9  qubits,  and  each  curve  extends  to  Tr/flflmd  for  the  corresponding  Ne. 


nHj)  =  I  I  pM(x,tMt)dt 
-Ml  J  0 


dx 


X^(t)dt 


dx. 


This  ratio  is  plotted  in  Fig.  2(b)  for  the  same  state  used 
above,  for  which  the  actual  ratio  is  one.  The  graph  shows 
that  reaching  the  actual  ratio  of  Ps  to  Pc  also  scales  exponen¬ 
tially  with  Ne.  However,  precise  values  for  measurement 
probabilities  are  often  not  needed  and  useful  qualitative  in¬ 
formation  can  be  obtained  for  integration  intervals  that  are 
shorter  than  7r/fi(und.  For  instance,  for  this  example  the 
curves  indicate  that  for  integration  limits  beyond  it/  (4fl(und), 
the  ratio  r  is  about  0.25,  of  the  same  order  as  the  actual 
value. 


V.  SIMULATION  EXAMPLE:  FACTORING 
A.  Quantum  Fourier  transform  and  Shor’s  algorithm 

The  most  celebrated  quantum  algorithm  is  Shor’s  method 
of  finding  the  prime  factors  of  a  number  N.  The  problem  of 
factorization  can  be  related  to  the  problem  of  determining  the 
period  p  of  the  function  f(x)=ax( mod  N),  for  an  integer  a 
that  is  co-prime  with  N;  if  p  is  even,  either  (ap,2+ 1)  or 
(apl2~  1)  will  have  a  common  factor  with  N  [8].  Shor’s  algo¬ 
rithm  relies  on  application  of  the  quantum  Fourier  transform 
(QFT)  to  f{x)  to  efficiently  determine  the  period. 

The  qubits  involved  in  implementing  this  algorithm  are 
divided  into  two  registers,  the  states  of  which  are  treated  as 
integers  according  to  the  states  of  the  associated  qubits.  The 
first  step  in  the  procedure  is  to  prepare  the  first  register  in  a 
superposition  of  all  computational  basis  states, 

2*1-1 

2  k), 

x=0 

by  applying  a  Hadamard  transform  to  each  of  the  N\  qubits 
in  the  register  [9].  A  gate  U  applied  to  both  registers  pro¬ 
duces  the  entangled  state 


2*1-1 

2  |x)  ®  |aA'(mod  N)).  (18) 

x=0 

If  the  second  register  is  measured,  it  collapses  to 
|aA°(mod  N))  for  some  x0.  The  first  register  ends  up  in  a 
superposition  of  all  states  |x')  for  which  ax  (mod  N) 
=  aA°(mod  N),  leaving  the  system  in  the  state 

(|xq)  +  |xq  +  p)  +  |x0  +  2 p)  +  ...)  ®  |av°(mod  N)).  (19) 

Application  of  the  QFT  to  the  first  register  imposes  interfer¬ 
ence  that  results  in  a  superposition  of  states  |x)  that  are  close 
to  integral  multiples  of  the  inverse  period,  x~xK=  k(2n'I p), 
for  integer  k.  Measurement  yields  an  integer  close  to  one  of 
the  xK,  and  after  several  iterations  the  period  p  can  be  deter¬ 
mined. 

The  minimum  number  of  qubits  required  for  each  register 
is  A' |  =  logo  N2  and  /V2=  log2  N.  This  ensures  that  the  second 
register  is  large  enough  to  represent  p,  which  satisfies  p 
<  N,  and  that  the  first  register  is  large  enough  to  give  a 
unique  value  for  p  from  the  QFT  (see  Ref.  [8]).  We  apply  our 
model  to  the  factorization  of  N  =  21  using  a  =  2,  the  first  in¬ 
teger  co-prime  with  21.  This  requires  14  qubits,  nine  for  the 
first  register  and  five  for  the  second;  the  algorithm  for  our 
example  is  illustrated  in  Fig.  3.  The  qubits  are  labeled  1 
through  14,  with  the  convention  that  qubit  1  (14)  is  the  most 
(least)  significant  qubit  for  the  first  (second)  register.  Qubit  n 
is  represented  using  the  frequency  a)n=(ol2n~1,  and  the  func¬ 
tion  \P  representing  the  state  of  the  system  is  defined  with  a 
resolution  of  1/w  over  an  interval  of  27r/(lfund=27r(2"max~1) 
[10].  In  Fig.  4  we  show  the  function  representing  the  state  of 
the  system  at  various  stages  in  the  calculation,  as  denoted  by 
the  dashed  lines  in  Fig.  3(a). 

We  implement  U  by  generating  the  output  state,  shown  in 
Fig.  4(b),  by  summing  all  functions  representing 
|x)|aA(mod  N))  over  x=0  to  x=214-l  [11].  From  here,  we 
measure  qubits  10  through  14  by  comparing  to 

pi")* .  F(n>  and  applying  the  rule  that  the  outcome  of  a  mea¬ 
surement  is  the  state  with  the  higher  measurement  probabil- 
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9 

10 
14 


FIG.  3.  (a)  Illustration  of  Shor’s  factoring  algorithm  for  N  =  21. 
The  first  gate  shown  represents  an  application  of  the  Hadamard 
transform  to  each  of  the  nine  qubits  in  the  first  register.  The  solid 
vertical  lines  that  terminate  the  circuit  flow  represent  measurement; 
the  dashed  vertical  lines  indicate  points  at  which  ’F  is  plotted  in  Fig. 
4.  (b)  Circuit  for  implementing  the  quantum  Fourier  transform 
shown  in  (a).  The  unitary  gate  Rd  is  a  rotation  by  the  angle  77/2  , 
and  H  is  a  Hadamard  transform.  The  QFT  reverses  the  order  of  the 
qubits. 


ity,  or  is  chosen  at  random  if  the  probabilities  are  equal.  This 
leaves  the  second  register  in  the  state  1 10000)=  1 16),  and  the 
first  register  in  a  linear  combination  of  all  |x')  between  0  and 
29-l  for  which  2X  (mod  21)=  16;  the  function  representing 
the  first  register  at  this  stage  is  shown  in  Fig.  4(c).  Qubits  10 
through  14  are  removed  from  the  calculation  and  the  QFT  is 
applied  to  the  remaining  nine  qubits. 

The  results  of  applying  the  QFT  and  subsequent  measure¬ 
ment  of  qubits  1-9  are  shown  in  Fig.  5(a),  which  displays  the 
measurement  probability  for  the  state  |x).  The  peaks  at  mul¬ 
tiples  of  85 y  indicate  a  period  of  p= 6  for  the  function 
2*(mod  21),  which  in  turn  yields  either  of  the  prime  factors 
of  21  from  gcd(ap/2+ 1 , N)  and  gcd(ap/2- 1 , N).  Figures 
5(b)-5(d)  zoom  in  on  the  probability  distributions  near  dif¬ 
ferent  multiples  of  85  j.  When  kX  85 y  is  not  an  integer,  the 
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FIG.  5.  (Color  online)  Results  of  the  QFT  applied  to  the  first 
register,  (a)  Plot  of  measurement  probability  for  all  x  between  0  and 
511.  The  six  peaks  are  designated  by  the  integer  k= 0  through  5. 
(b)-(d)  Details  of  peaks  for  k=  1,  2  and  3  in  (a).  The  dashed  line 
indicates  the  exact  value  of  k(2NiIp)  for  that  peak,  (e)  Plot  of 
measurement  probability  versus  jc'  for  the  case  when  the  QFT  is  not 
used.  Only  part  of  the  entire  distribution,  which  extends  to  x' 
=  511,  is  shown. 

probability  is  distributed  among  integers  closest  to  the  frac¬ 
tional  value,  as  in  (b)  and  (c). 

B.  Simplifications  in  a  classical  model 

There  is  a  dramatic  simplification  to  the  factoring  algo¬ 
rithm  when  using  a  classical  model  for  simulation.  The  need 
for  the  QFT  stems  from  the  fact  that,  in  a  quantum  system, 
measurement  of  the  state  in  Eq.  (19)  results  in  one  of  the 
states  in  the  superposition,  with  no  opportunity  to  learn  the 
others.  Repetition  of  the  algorithm  to  that  point  likely  results 
in  a  different  x0,  preventing  the  period  from  being  learned  in 
successive  iterations.  In  a  classical  system,  the  state  in  Eq. 
(19)  can  be  measured  as  many  times  as  necessary  to  deter¬ 
mine  the  period  p.  If  we  apply  this  simplified  scheme  to  the 
state  in  Fig.  4(c),  we  find  equal  probability  for  any  x' 
satisfying  2X  (mod  21)=  16;  the  first  eight  peaks  are  shown  in 
Fig.  5(e). 

In  addition  to  avoiding  all  of  the  gates  required  for  the 
QFT,  no  imaginary  numbers  are  required,  and  most  impor¬ 
tantly,  we  save  on  qubits.  The  first  register,  whose  size  for 
Shor’s  algorithm  is  dictated  by  the  QFT  stage,  in  this  case 
only  needs  enough  qubits  to  represent  p,  and  p<N.  This 
provides  a  savings  on  order  of  log2  N  qubits  compared  to  the 
quantum  case.  Applying  this  simplified  factoring  algorithm 
to  N=  21  would  require  only  10  qubits,  five  for  each  register. 


FIG.  4.  (Color  online)  Function  'F  representing  the  state  of  the 
system  at  various  stages  in  the  simulation  plotted  over  one  period. 
The  different  durations  shown  are  due  to  the  different  number  of 
qubits  at  those  points  in  the  calculation.  Only  the  real  part  of  the 
function  is  plotted  in  (d). 


VI.  CONCLUSION 

We  have  presented  a  classical,  qubit-based  model  of  dis¬ 
crete  quantum  systems  that  offers  a  framework  for  simulat¬ 
ing  quantum  computations  by  providing  access  to  individual 
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qubits  and  their  interactions.  The  dimension  of  the  state 
space  for  a  composite  system  grows  exponentially  with  qubit 
number,  and  individual  qubits  continue  to  be  accessible.  Ap¬ 
plication  to  Shor’s  algorithm  highlights  the  features  of  the 
model  in  action.  The  resources  required  for  implementation 
scale  exponentially  with  the  number  of  entangled  qubits,  yet 
it  is  possible  to  save  on  qubits  in  a  classical  model. 
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APPENDIX:  FOURIER  DECOMPOSITION  OF  HN(t) 

The  HN{t)  can  be  Fourier  decomposed  by  expanding  the 
products  of  harmonics  hn  in  terms  of  sum  and  difference 
frequencies.  Here  we  explicitly  show  the  decomposition  for 
the  case  of  three  qubits.  The  eight  functions  H2j=1  &(t)  can 
be  expanded  as 


cos(w1f)cos(w20cos(w3?)  =  “(cos[(&q  +  <w2  +  «3)f]  +  cos[(«i  +  a>2  -  «3)f]  +  cos[(wj  -  a>2  +  w3)r]  +  cos[(aq  -  a>2  -  cu3)f]), 

cos(w1f)cos(«w20sin(<i)3f)  =  “(sin[(wi  +  u>2+  &>3)f]  -  sin[(&q  +  u>2  -  cu3)f]  +  sin[(oq  -  w2  +  co3)f]  -  sin[(w1  -  a>2—  o>3)f]). 


cos(o)1f)sin(a)2f)cos(w3f)  =  -(sin[(wj  +  co2  +  «3)f]  +  sin[(aq  +  u>2  -  <w3)f]  -  sinfCcO]  -  u>2  +  <w3)f]  -  sin[(&q  -  co2  -  co3)f]), 


cos  ( u) !  t)  sin(  cj2t)  sin(  w3 1) 


1 

4 


(-  cos[(w1  +  a>2  +  w3)f]  +  cos[(o)|  +  a>2  -  w3)r]  +  cos[(oq  - 


Ct)2  +  <W3)f]  -  COs[(«]  -  co2—  w3)t]). 


sin(<u1f)cos(w2f)cos (co3f)  =  -(sin[(co1  +  co2  +  «3)f]  +  sin[(w1  +  u>2  -  cu3)f]  +  sin[(a>1  -  w2  +  w3)f]  +  sin[(tiq  -  co2  -  cn3)f]), 


sin(w1f)cos(<u2f)sin(<u3f)  =  -(-  cos[(&q  +  co2  +  w3)f]  +  cos^oq  +  a>2-  <w3)f]  -  cos[(oq  -  co2  +  at 3)f]  +  cos[(oq  -  w2  -  o>3)t]), 

sin(w1f)sin(w2f)cos((u3f)  =  -(-  cos[(oq  +  co2  +  w3)f]  -  cos[(oq  +  a>2  -  <w3)r]  +  cos[(oq  -  co2  +  &>3)r]  +  cos[(oq  —  (o2—  &>3)f]), 

sin(w1f)sin(w2f)sin(w3?)  =  —  (—  sin[(oq  +  u>2  +  &)3)f]  +  sin[(oq  +  <w2  -  &>3)f]  +  sin[(aq  -  a>2  +  «u3)f]  -  sin[(oq  -  ut2  -  w3)f]) . 


(Al) 


We  introduce  a  more  compact  notation  to  represent  the  eight 
Fourier  components: 

(1,0,0, 0,0,0, 0,0)  =  -cos[(oq  +  co2  +  «3)f] 

(0, 1,0, 0,0,0, 0,0)  =  -cos[(oq  +  (n)2-  w3)f] 

(0,0, 1,0, 0,0, 0,0)  =  -cos[(oq  -  co2  +  «3)f] 

(0,0,0, 1,0,0, 0,0)  =  -cos[(oq  -  co2  -  w3)f] 


(0, 0,0,0, 1,0, 0,0)  =  -sin[(oq  +  w2  +  <w3)f] 

(0,0,0, 0,0, 1,0,0)  =  -sin[(oq  +  w2  -  w3)f] 

(0,0,0, 0,0,0, 1,0)  =  -sin[(oq  -  w2  +  <w3)f] 

(0,0, 0,0,0, 0,0, 1)  =  -sin[(uq  -  w2  -  «3)f].  (A2) 

For  all  of  these  terms,  ( 1  •  ( 1  tt/  1 60fund) ,  where  (lj 

corresponds  to  the  Fourier  component  represented  by  the 
row  vector  with  a  one  in  the  ath  place  and  zeros  everywhere 
else. 
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The  functions  //3(f)  written  with  this  vector  notation  be¬ 
come 

cos(«if)cos(ci>2f)cos(w3f)  =  (1,1, 1,1, 0,0, 0,0), 

cos(a)1f)cos(w2f)sin(w3f)  =  (0,0, 0,0, 1,-  1,1,—  1), 
cos(w1f)sin(a)2f)cos(ai3f)  =  (0,0, 0,0, 1, 1,-  1,—  1), 
cos(<u1f)sin(<u2f)sin((U3f)  =  (-  1,1,1,-  1,0,0, 0,0), 
sin(<w1f)cos(<w2f)cos((U3f)  =  (0,0, 0,0, 1, 1, 1, 1), 
sin(<xi1f)cos(&)2f)sin(&)3f)  =  (-  1,1,-  1, 1,0,0, 0,0), 
sin(<x>1f)sin(<y2f)cos(«3f)  =  (-  1,—  1, 1, 1,0,0, 0,0), 

sin(w1f)sin(o)2f)sin(o)3f)  =  (0,0, 0,0,-  1,1,1,—  1). 

From  this  it  can  be  seen  that  all  of  the  H2{t)  are  orthogonal. 
For  example. 


cos(<w1f)cos(<w2f)cos(<w3f)  •  cos(w1f)sin(w2f)sin(w3f) 

=  (1,1,1,1,0,0,0,0)  •  (-  1,1,1,-  1,0,0,0,0) 

=  -1  +  1  +  1- 1=0.  (A3) 

In  this  case  and  for  small  numbers  of  qubits,  all  of  the  dif¬ 
ferent  inner  products  can  be  verified  to  be  zero,  and 
Hnj  '  ^«j=(7r/4Iifund). 

We  can  also  see  that  not  all  of  the  functions  are  orthogo¬ 
nal  if  there  is  redundancy  in  the  Fourier  frequencies.  If,  for 
instance,  u>l  +  a)2-  co2  =  col  -  u>2+  u>2,  then  the  second  and  third 
vectors  in  Eq.  (A2)  are  identical  (as  are  the  sixth  and  sev¬ 
enth).  In  this  case,  the  H2{t)  can  only  be  represented  by  six 
orthogonal  vectors,  not  eight.  Then, 

cos(&>1f)cos(<w2f)cos(<w3f)  •  cos(<Uif)sin(a>2f)sin(<w3f) 

=  (1, 2, 1, 0,0,0) -(-1, 2,-1, 0,0,0) 

=  -  1  +  4  -  1  =  2.  (A4) 
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